#include "phylip.h"

/* version 3.56c. (c) Copyright 1993 by Joseph Felsenstein.
   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
   Permission is granted to copy and use this program provided no fee is
   charged for it and provided that this copyright notice is not removed. */

#define namelength      10   /* number of characters in species name    */
#define epsilon         0.000001   /* a very small but not too small number */

#define ibmpc0          false
#define ansi0           true
#define vt520           false
#define down            2
#define over            60

typedef double *vector;		  /* nodes will form binary tree           */
typedef char naym[namelength];

struct node
{		  /* describes a tip species or an ancestor */
	node*	next;
	node*	back;
	long	index;
	bool	tip;				/* present species are tips of tree       */
	vector	d;					/* distances							  */
	vector	w;					/* weights				                  */
	double	t;					/* time                                   */
	bool	sametime;			/* bookkeeps scrunched nodes              */
	double	weight;				/* weight of node used by scrunch         */
	bool	processed;			/* used by evaluate                       */
};

typedef node **pointptr;
typedef long longer[6];

node *root, *best;
FILE *infile, *outfile, *treefile;
long numsp, numsp2, inseed, col, datasets, ith,
i, j, l, jumb, njumble;
/* numsp = number of species*/
pointptr treenode, bestree;	  /* pointers to all nodes in tree */
vector<string> nayms;
bool jumble, lower, upper, replicates, trout,
printdata, treeprint, mulsets, ibmpc, vt52,
ansi, firstset;
double power;
longer seed;
long *enterorder;
char ch;
/* Local variables for maketree, propagated globally for C version: */
long examined;
double like, bestyet;
node *there;


void getnums()
{
	/* read species number */
	fscanf(infile, "%ld", &numsp);
	fprintf(outfile, "\n%4ld Populations\n", numsp);
	numsp2 = numsp * 2 - 1;
}  /* getnums */

void getoptions()
{
	/* interactively set options */
	long i, inseed0=0;
	char ch;
	bool done=false, done1=false;

	fprintf(outfile, "\nFitch-Margoliash method ");
	fprintf(outfile, "with contemporary tips, version %s\n\n",VERSION);
	jumble = false;
	njumble = 1;
	lower = false;
	power = 2.0;
	replicates = false;
	upper = false;
	trout = true;
	printdata = false;
	treeprint = true;
	for (;;)
	{
		printf( ansi ? "\033[2J\033[H" :
				vt52 ? "\033E\033H"    :
				"\n");
		printf("\nFitch-Margoliash method ");
		printf("with contemporary tips, version %s\n\n",VERSION);
		printf("Settings for this run:\n");
		printf("  U                 Search for best tree?  %s\n",
			   usertree ? "No, use user trees in input file" : "Yes");
		printf("  P                                Power?%9.5f\n",power);
		printf("  L         Lower-triangular data matrix?  %s\n",
			   (lower ? "Yes" : "No"));
		printf("  R         Upper-triangular data matrix?  %s\n",
			   (upper ? "Yes" : "No"));
		printf("  S                        Subreplicates?  %s\n",
			   (replicates ? "Yes" : "No"));
//		if (!usertree)
		{
			printf("  J     Randomize input order of species?");
			if (jumble)
				printf("  Yes (seed =%8ld,%3ld times)\n", inseed0, njumble);
			else
				printf("  No. Use input order\n");
		}
		printf("  M           Analyze multiple data sets?");
		if (mulsets)
			printf("  Yes, %2ld sets\n", datasets);
		else
			printf("  No\n");
		printf("  0   Terminal type (IBM PC, VT52, ANSI)?  %s\n",
			   (ibmpc ? "IBM PC" :
				ansi  ? "ANSI"   :
				vt52  ? "VT52"   :
				"(none)"));

		printf("  1    Print out the data at start of run  %s\n",
			   (printdata ? "Yes" : "No"));
		printf("  3                        Print out tree  %s\n",
			   (treeprint ? "Yes" : "No"));
		printf("  4       Write out trees onto tree file?  %s\n",
			   (trout ? "Yes" : "No"));
		printf("\nAre these settings correct?");
		printf(" (type Y or the letter for one to change)\n");
		scanf("%c%*[^\n]", &ch);
		getchar();
		if (ch == '\n')
			ch = ' ';
		uppercase(&ch);
		if (ch == 'Y')
			break;
		if (strchr("JUP-LRSM12340",ch))
		{
			switch (ch)
			{
				
				case 'J':
					jumble = !jumble;
					if (jumble)
					{
						printf("Random number seed (must be odd)?\n");
						scanf("%ld%*[^\n]", &inseed);
						getchar();
						inseed0 = inseed;
						for (i = 0; i <= 5; i++)
							seed[i] = 0;
						i = 0;
						do
						{
							seed[i] = inseed & 63;
							inseed /= 64;
							i++;
						} while (inseed != 0);
						printf("Number of times to jumble?\n");
						scanf("%ld%*[^\n]", &njumble);
						getchar();
					}
					else njumble = 1;
					break;

				case 'L':
					lower = !lower;
					break;

				case 'P':
					printf("New power?\n");
					scanf("%lf%*[^\n]", &power);
					getchar();
					break;

				case 'R':
					upper = !upper;
					break;

				case 'S':
					replicates = !replicates;
					break;


				case 'M':
					mulsets = !mulsets;
					if (mulsets)
					{
						done1 = false;
						do
						{
							printf("How many data sets?\n");
							scanf("%ld%*[^\n]", &datasets);
							getchar();
							done1 = (datasets >= 1);
							if (!done1)
								printf("BAD DATA SETS NUMBER:  it must be greater than 1\n");
						} while (done1 != true);
					}
					break;

				case '0':
					if (ibmpc)
					{
						ibmpc = false;
						vt52 = true;
					}
					else
					{
						if (vt52)
						{
							vt52 = false;
							ansi = true;
						}
						else if (ansi)
							ansi = false;
						else
							ibmpc = true;
					}
					break;

				case '1':
					printdata = !printdata;
					break;

				case '3':
					treeprint = !treeprint;
					break;

				case '4':
					trout = !trout;
					break;
			}
		}
		else
			printf("Not a possible option!\n");
	}
}  /* getoptions */


void doinit()
{
	/* initializes variables */
	long i, j;
	node *p, *q;

	getnums();
	getoptions();
	treenode = (node **)Malloc(numsp2*sizeof(node *));
	for (i = 0; i < (numsp); i++)
	{
		treenode[i] = (node *)Malloc(sizeof(node));
		treenode[i]->d = (vector)Malloc(numsp2*(sizeof(double)));
		treenode[i]->w = (vector)Malloc(numsp2*(sizeof(double)));
	}
	for (i = numsp; i < (numsp2); i++)
	{
		q = NULL;
		for (j = 1; j <= 3; j++)
		{
			p = (node *)Malloc(sizeof(node));
			p->d = (vector)Malloc(numsp2*(sizeof(double)));
			p->w = (vector)Malloc(numsp2*(sizeof(double)));
			p->next = q;
			q = p;
		}
		p->next->next->next = p;
		treenode[i] = p;
	}
//	if (!usertree && njumble > 1)
	bestree = (node **)Malloc(numsp2*sizeof(node *));
	for (i = 0; i < (numsp); i++)
	{
		bestree[i] = (node *)Malloc(sizeof(node));
		bestree[i]->d = (vector)Malloc(numsp2*(sizeof(double)));
		bestree[i]->w = (vector)Malloc(numsp2*(sizeof(double)));
	}
	for (i = numsp; i < (numsp2); i++)
	{
		q = NULL;
		for (j = 1; j <= 3; j++)
		{
			p = (node *)Malloc(sizeof(node));
			p->d = (vector)Malloc(numsp2*(sizeof(double)));
			p->w = (vector)Malloc(numsp2*(sizeof(double)));
			p->next = q;
			q = p;
		}
		p->next->next->next = p;
		bestree[i] = p;
	}

}  /* doinit */

void inputoptions()
{
	/* read options information */
	char ch;
	long cursp;

	if (!firstset)
	{
		if (eoln(infile))
		{
			fscanf(infile, "%*[^\n]");
			getc(infile);
		}
		fscanf(infile, "%ld", &cursp);
		if (cursp != numsp)
		{
			printf("\nERROR: INCONSISTENT NUMBER OF SPECIES IN DATA SET %4ld\n",
				   ith);
			exit(-1);
		}
	}
	while (!eoln(infile))
	{
		ch = getc(infile);
		if (ch == '\n')
			ch = ' ';
		uppercase(&ch);
		if (ch != ' ')
		{
			printf("BAD OPTION CHARACTER: %c\n", ch);
			exit(-1);
		}
	}
	fprintf(outfile, "                  __ __             2\n");
	fprintf(outfile, "                  \\  \\   (Obs - Exp)\n");
	fprintf(outfile, "Sum of squares =  /_ /_  ------------\n");
	fprintf(outfile, "                               ");
	if (power == (long)power)
		fprintf(outfile, "%2ld\n", (long)power);
	else
		fprintf(outfile, "%4.1f\n", power);
	fprintf(outfile, "                   i  j      Obs\n\n");
}  /* inputoptions */



void getdata()
{
	/* read in distance matrix */
	long i, j, k, columns, n;
	bool skipit, skipother;
	double x;
	node *p;

	columns = replicates ? 4 : 6;
	if (printdata)
	{
		fprintf(outfile, "\nName                       Distances");
		if (replicates)
			fprintf(outfile, " (replicates)");
		fprintf(outfile, "\n----                       ---------");
		if (replicates)
			fprintf(outfile, "-------------");
		fprintf(outfile, "\n\n");
	}
	for (i = 1; i <= (numsp2); i++)
	{
		treenode[i - 1]->back = NULL;
		treenode[i - 1]->index = i;
		treenode[i - 1]->tip = (i <= numsp);
		treenode[i - 1]->t = 0.0;
		treenode[i - 1]->sametime = false;
		if (i > numsp)
		{
			p = treenode[i - 1]->next;
			while (p != treenode[i - 1])
			{
				p->back = NULL;
				p->tip = false;
				p->index = i;
				p = p->next;
			}
		}
	}
	if (njumble > 1)
		for (i = 1; i <= (numsp2); i++)
		{
			bestree[i - 1]->back = NULL;
			bestree[i - 1]->index = i;
			bestree[i - 1]->tip = (i <= numsp);
			bestree[i - 1]->t = 0.0;
			bestree[i - 1]->sametime = false;
			if (i > numsp)
			{
				p = bestree[i - 1]->next;
				while (p != bestree[i - 1])
				{
					p->back = NULL;
					p->tip = false;
					p->index = i;
					p = p->next;
				}
			}
		}
	for (i = 0; i < (numsp); i++)
	{
		treenode[i]->d[i] = 0.0;
		treenode[i]->w[i] = 0.0;
		treenode[i]->weight = 0.0;
		infile >> names[i];
		for (j = 1; j <= (numsp); j++)
		{
			skipit = ((lower && j >= i + 1) || (upper && j <= i + 1));
			skipother = ((lower && i + 1 >= j) || (upper && i + 1 <= j));
			if (!skipit)
			{
				if (eoln(infile))
				{
					fscanf(infile, "%*[^\n]");
					getc(infile);
				}
				fscanf(infile, "%lf", &x);
				treenode[i]->d[j - 1] = x;
				if (replicates)
				{
					if (eoln(infile))
					{
						fscanf(infile, "%*[^\n]");
						getc(infile);
					}
					fscanf(infile, "%ld", &n);
				}
				else
					n = 1;
				if (n > 0 && x < 0)
				{
					printf("NEGATIVE DISTANCE BETWEEN SPECIES%5ld AND %5ld\n",
						   i + 1, j);
					exit(-1);
				}
				treenode[i]->w[j - 1] = n;
				if (skipother)
				{
					treenode[j - 1]->d[i] = treenode[i]->d[j - 1];
					treenode[j - 1]->w[i] = treenode[i]->w[j - 1];
				}
			}
		}
	}
	for (i = 0; i < (numsp); i++)
	{
		for (j = 0; j < (numsp); j++)
		{
			if (i + 1 != j + 1)
			{
				if (treenode[i]->d[j] < epsilon)
					treenode[i]->d[j] = epsilon;
				treenode[i]->w[j] /= exp(power * log(treenode[i]->d[j]));
			}
		}
	}
}  /* getdata */

void add(node * below, node *newtip, node * newfork)
{
	/* inserts the nodes newfork and its left descendant, newtip,
	   to the tree.  below becomes newfork's right descendant */
	if (below != treenode[below->index - 1])
		below = treenode[below->index - 1];
	if (below->back != NULL)
		below->back->back = newfork;
	newfork->back = below->back;
	below->back = newfork->next->next;
	newfork->next->next->back = below;
	newfork->next->back = newtip;
	newtip->back = newfork->next;
	if (root == below)
		root = newfork;
	root->back = NULL;
}  /* add */

void re_move(node **item, node **fork)
{
	/* removes nodes item and its ancestor, fork, from the tree.
	   the new descendant of fork's ancestor is made to be
	   fork's second descendant (other than item).  Also
	   returns pointers to the deleted nodes, item and fork */
	node *p, *q;

	if ((*item)->back == NULL)
	{
		*fork = NULL;
		return;
	}
	*fork = treenode[(*item)->back->index - 1];
	if (root == *fork)
	{
		if (*item == (*fork)->next->back)
			root = (*fork)->next->next->back;
		else
			root = (*fork)->next->back;
	}
	p = (*item)->back->next->back;
	q = (*item)->back->next->next->back;
	if (p != NULL)
		p->back = q;
	if (q != NULL)
		q->back = p;
	(*fork)->back = NULL;
	p = (*fork)->next;
	while (p != *fork)
	{
		p->back = NULL;
		p = p->next;
	}
	(*item)->back = NULL;
}  /* remove */

void scrunchtraverse(node *u, node  **closest, double *tmax)
{
	/* traverse to find closest node to the current one */
	if (!u->sametime)
	{
		if (u->t > *tmax)
		{
			*closest = u;
			*tmax = u->t;
		}
		return;
	}
	u->t = treenode[u->back->index - 1]->t;
	if (!u->tip)
	{
		scrunchtraverse(u->next->back, closest,tmax);
		scrunchtraverse(u->next->next->back, closest,tmax);
	}
}  /* scrunchtraverse */

void combine(node *a, node *b)
{
	/* put node b into the set having the same time as a */
	if (a->weight + b->weight <= 0.0)
		a->t = 0.0;
	else
		a->t = (a->t * a->weight + b->t * b->weight) / (a->weight + b->weight);
	a->weight += b->weight;
	b->sametime = true;
}  /* combine */

void scrunch(node *s)
{
	/* see if nodes can be combined to prevent negative lengths */
/* Local variables for scrunch: */
	
	bool found;
	node *closest = 0;
	double tmax = -1.0;
	do
	{
		if (!s->tip)
		{
			scrunchtraverse(s->next->back, &closest,&tmax);
			scrunchtraverse(s->next->next->back, &closest,&tmax);
		}
		found = (tmax > s->t);
		if (found)
			combine(s, closest);
		tmax = -1.0;
	} while (found);
}  /* scrunch */

void secondtraverse(node *a, node *q, node* u, node * v, long i, long j, long k, double *sum)
{
	/* recalculate distances, add to sum */
	long l;
	double wil, wjl, wkl, wli, wlj, wlk, TEMP;

	if (!(a->processed || a->tip))
	{
		secondtraverse(a->next->back, q,u,v,i,j,k,sum);
		secondtraverse(a->next->next->back, q,u,v,i,j,k,sum);
		return;
	}
	if (!(a != q && a->processed))
		return;
	l = a->index;
	wil = u->w[l - 1];
	wjl = v->w[l - 1];
	wkl = wil + wjl;
	wli = a->w[i - 1];
	wlj = a->w[j - 1];
	wlk = wli + wlj;
	q->w[l - 1] = wkl;
	a->w[k - 1] = wlk;
	if (wkl <= 0.0)
		q->d[l - 1] = 0.0;
	else
		q->d[l - 1]	= (wil * u->d[l - 1] + wjl * v->d[l - 1]) / wkl;
	if (wlk <= 0.0)
		a->d[k - 1] = 0.0;
	else
		a->d[k - 1]	= (wli * a->d[i - 1] + wlj * a->d[j - 1]) / wlk;
	if (wkl > 0.0)
	{
		TEMP = u->d[l - 1] - v->d[l - 1];
		(*sum) += wil * wjl / wkl * (TEMP * TEMP);
	}
	if (wlk > 0.0)
	{
		TEMP = a->d[i - 1] - a->d[j - 1];
		(*sum) += wli * wlj / wlk * (TEMP * TEMP);
	}
}  /* secondtraverse */

void firstraverse(node *q_, node *r, double *sum)
{  /* firsttraverse                              */
   /* go through tree calculating branch lengths */
	/* Local variables for firstraverse: */
	node *q;
	long i, j, k;
	node *u, *v;

	q = q_;
	if (q == NULL)
		return;
	q->sametime = false;
	if (!q->tip)
	{
		firstraverse(q->next->back, r,sum);
		firstraverse(q->next->next->back, r,sum);
	}
	q->processed = true;
	if (q->tip)
		return;
	u = q->next->back;
	v = q->next->next->back;
	i = u->index;
	j = v->index;
	k = q->index;
	if (u->w[j - 1] + v->w[i - 1] <= 0.0)
		q->t = 0.0;
	else
		q->t = (u->w[j - 1] * u->d[j - 1] +
				v->w[i - 1] * v->d[i - 1]) /
			   (2.0 * (u->w[j - 1] + v->w[i - 1]));
	q->weight = u->weight + v->weight + u->w[j - 1] + v->w[i - 1];
	scrunch(q);
	secondtraverse(r,q,u,v,i,j,k,sum);
}  /* firstraverse */

void sumtraverse(node *q, double *sum)
{
	/* traverse to finish computation of sum of squares */
	long i, j;
	node *u, *v;
	double TEMP, TEMP1;

	if (q->tip)
		return;
	sumtraverse(q->next->back, sum);
	sumtraverse(q->next->next->back, sum);
	u = q->next->back;
	v = q->next->next->back;
	i = u->index;
	j = v->index;
	TEMP = u->d[j - 1] - 2.0 * q->t;
	TEMP1 = v->d[i - 1] - 2.0 * q->t;
	(*sum) += u->w[j - 1] * (TEMP * TEMP) + v->w[i - 1] * (TEMP1 * TEMP1);
}  /* sumtraverse */

void evaluate(node *r)
{
	/* fill in times and evaluate sum of squares for tree */
	/* Local variables for evaluate: */
	double sum;
	long i;
	sum = 0.0;
	for (i = 0; i < (numsp2); i++)
		treenode[i]->processed = treenode[i]->tip;
	firstraverse(r, r,&sum);
	sumtraverse(r, &sum);
	examined++;
	if (replicates && (lower || upper))
		sum /= 2;
	like = -sum;
}  /* evaluate */


void tryadd(node *p, node **item, node **nufork)
{
	/* temporarily adds one fork and one tip to the tree.
	   if the location where they are added yields greater
	   "likelihood" than other locations tested up to that
	   time, then keeps that location as there */
	add(p, *item, *nufork);
	evaluate(root);
	if (like > bestyet)
	{
		bestyet = like;
		there = p;
	}
	re_move(item, nufork);
}  /* tryadd */

void addpreorder(node *p, node *item, node *nufork)
{
	/* traverses a binary tree, calling PROCEDURE tryadd
	   at a node before calling tryadd at its descendants */
/* Local variables for addpreorder: */
	if (p == NULL)
		return;
	tryadd(p, &item, &nufork);
	if (!p->tip)
	{
		addpreorder(p->next->back, item, nufork);
		addpreorder(p->next->next->back, item, nufork);
	}
}  /* addpreorder */

void tryrearr(node *p, node **r, bool *success)
{
	/* evaluates one rearrangement of the tree.
	   if the new tree has greater "likelihood" than the old
	   one sets success := TRUE and keeps the new tree.
	   otherwise, restores the old tree */
	node *frombelow, *whereto, *forknode;
	double oldlike;

	if (p->back == NULL)
		return;
	forknode = treenode[p->back->index - 1];
	if (forknode->back == NULL)
		return;
	oldlike = like;
	if (p->back->next->next == forknode)
		frombelow = forknode->next->next->back;
	else
		frombelow = forknode->next->back;
	whereto = forknode->back;
	re_move(&p, &forknode);
	add(whereto, p, forknode);
	if ((*r)->back != NULL)
		*r = treenode[(*r)->back->index - 1];
	evaluate(*r);
	if (like > oldlike)
	{
		bestyet = like;
		*success = true;
		return;
	}
	re_move(&p, &forknode);
	add(frombelow, p, forknode);
	if ((*r)->back != NULL)
		*r = treenode[(*r)->back->index - 1];
	like = oldlike;
}  /* tryrearr */

void repreorder(node *p, node **r, bool *success)
{
	/* traverses a binary tree, calling PROCEDURE tryrearr
	   at a node before calling tryrearr at its descendants */
	if (p == NULL)
		return;
	tryrearr(p,r,success);
	if (!p->tip)
	{
		repreorder(p->next->back,r,success);
		repreorder(p->next->next->back,r,success);
	}
}  /* repreorder */

void rearrange(node **r_)
{
	/* traverses the tree (preorder), finding any local
	   rearrangement which decreases the number of steps.
	   if traversal succeeds in increasing the tree's
	   "likelihood", PROCEDURE rearrange runs traversal again */
/* Local variables for rearrange: */
	node **r;
	bool success;
	r = r_;
	success = true;
	while (success)
	{
		success = false;
		repreorder(*r,r,&success);
	}
}  /* rearrange */


void treeout(node *p)
{
	/* write out file with representation of final tree */
	long i, n, w;
	char c;
	double x;

	if (p->tip)
		treefile << nayms[p->index - 1];
	else
	{
		putc('(', treefile);
		col++;
		treeout(p->next->back);
		putc(',', treefile);
		col++;
		if (col > 55)
		{
			putc('\n', treefile);
			col = 0;
		}
		treeout(p->next->next->back);
		putc(')', treefile);
		col++;
	}
	if (p != root)
		x = treenode[p->back->index - 1]->t - p->t;
	if (x > 0.0)
		w = (long)(0.43429448222 * log(x));
	else if (x == 0.0)
		w = 0;
	else
		w = (long)(0.43429448222 * log(-x)) + 1;
	if (w < 0)
		w = 0;
	if (p == root)
		fprintf(treefile, ";\n");
	else
	{
		fprintf(treefile, ":%*.5f", (int)(w + 7), x);
		col += w + 8;
	}
}  /* treeout */


void copynode(node *c, node *d)
{
	/* make a copy of a node */

	memcpy(d->d, c->d, numsp2*sizeof(double));
	memcpy(d->w, c->w, numsp2*sizeof(double));
	d->t = c->t;
	d->sametime = c->sametime;
	d->weight = c->weight;
	d->processed = c->processed;
}  /* copynode */

void copy_(pointptr a, pointptr b)
{
	/* make a copy of a tree */
	short i, j=0;
	node *p, *q;

	for (i = 0; i < numsp; i++)
	{
		copynode(a[i], b[i]);
		if (a[i]->back != NULL)
		{
			if (a[i]->back == a[a[i]->back->index - 1])
				b[i]->back = b[a[i]->back->index - 1];
			else if (a[i]->back == a[a[i]->back->index - 1]->next)
				b[i]->back = b[a[i]->back->index - 1]->next;
			else
				b[i]->back = b[a[i]->back->index - 1]->next->next;
		}
		else b[i]->back	= NULL;
	}
	for (i = numsp; i < numsp2; i++)
	{
		p = a[i];
		q = b[i];
		for (j = 1; j <= 3; j++)
		{
			copynode(p, q);
			if (p->back)
			{
				if (p->back == a[p->back->index - 1])
					q->back = b[p->back->index - 1];
				else if (p->back == a[p->back->index - 1]->next)
					q->back = b[p->back->index - 1]->next;
				else
					q->back	= b[p->back->index - 1]->next->next;
			}
			else
				q->back	= NULL;
			p = p->next;
			q = q->next;
		}
	}
}  /* copy */


void maketree()
{
	/* constructs a binary tree from the pointers in treenode.
	   adds each node at location which yields highest "likelihood"
	   then rearranges the tree for greatest "likelihood" */
	long i, j, k;
	double bestlike, bstlike2, gotlike;
	bool lastrearr;
	node *item, *nufork;

	if (jumb == 1)
	{
		getdata();
		examined = 0;
	}
	for (i = 0; i < (numsp); i++)
		enterorder[i] = i;
	if (jumble)
	{
		for (i = 0; i < (numsp); i++)
			swap(enterorder[rand() % numsp], enterorder[i]);
	}
	root = treenode[enterorder[0]];
	add(treenode[enterorder[0]], treenode[enterorder[1]],
		treenode[numsp]);
	lastrearr = false;
	for (i = 3; i <= (numsp); i++)
	{
		bestyet = -10000000.0;
		item = treenode[enterorder[i - 1]];		
		nufork = treenode[numsp + i - 2];
		addpreorder(root, item, nufork);
		add(there, item, nufork);
		like = bestyet;
		rearrange(&root);
		evaluate(root);
		examined--;
		lastrearr = (i == numsp);
		if (lastrearr)
		{
			bestlike = bestyet;
			do
			{
				gotlike = bestlike;
				for (j = 0; j < (numsp2); j++)
				{
					there = root;
					bestyet = -32000.0;
					item = treenode[j];
					if (item != root)
					{
						re_move(&item, &nufork);
						there = root;
						addpreorder(root, item, nufork);
						add(there, item, nufork);
					}
				}
			} while (bestlike > gotlike);
			if (njumble > 1)
			{
				if (jumb == 1 || (jumb > 1 && bestlike > bstlike2))
				{
					copy_(treenode, bestree);
					best = bestree[root->index -1];
					bstlike2 = bestlike;
				}
			}
		}
		if (i == numsp && njumble == jumb)
		{
			if (njumble > 1)
			{
				copy_(bestree, treenode);
				root = treenode[best->index - 1];
			}
			evaluate(root);
//			printree();
//			describe();
		}
	}
}  /* maketree */


int main(int argc,  char *argv[])
{  /* Fitch-Margoliash criterion with contemporary tips */
	char infilename[100],outfilename[100],trfilename[100];
	/* reads in numsp, options, and the data, then calls maketree to
	   construct the tree */
	openfile(&infile,INFILE,"r",argv[0],infilename);
	openfile(&outfile,OUTFILE,"w",argv[0],outfilename);

	ibmpc = ibmpc0;
	ansi = ansi0;
	vt52 = vt520;
	mulsets = false;
	firstset = true;
	datasets = 1;
	doinit();
	openfile(&treefile,TREEFILE,"w",argv[0],trfilename);
	enterorder = (long *)Malloc(numsp*sizeof(long));
	for (ith = 1; ith <= datasets; ith++)
	{
		if (datasets > 1)
		{
			fprintf(outfile, "\nData set # %ld:\n",ith);
		}
		
		/* reads the input data */
		inputoptions();
		
		for (jumb = 1; jumb <= njumble; jumb++)
			maketree();
		firstset = false;
		if (eoln(infile))
		{
			fscanf(infile, "%*[^\n]");
			getc(infile);
		}
	}
	FClose(infile);
	FClose(outfile);
	FClose(treefile);
	exit(0);
}  /* Fitch-Margoliash criterion with contemporary tips */

